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We study the phase diagram of a generahzed Winfree model. The modification is such that 
the coupling depends on the fraction of synchronized oscillators, a situation which has been noted 
in some experiments on coupled Josephson junctions and mechanical systems. We let the global 
OO ■ coupling fc be a function of the Kuramoto order parameter r through an exponent z such that z = 1 

^^ ' corresponds to the standard Winfree model, z < 1 strengthens the coupling at low r (low amount 

^^ \ of synchronization) and, at 2 > 1, the coupling is weakened for low r. Using both analytical and 

Cn . numerical approaches, we find that 2 controls the size of the incoherent phase region, and one may 

make the incoherent behavior less typical by choosing 2 < 1. We also find that the original Winfree 
model is a rather special case, indeed the partial locked behavior disappears for 2 > 1. At fixed 
k and varying 7, the stability boundary of the locked phase corresponds to a transition that is 
^^ ' continuous for 2 < 1 and first-order for 2 > 1. This change in the nature of the transition is in 

Cm ' accordance with a previous study on a similarly modified Kuramoto model. 
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I. INTRODUCTION 



Forty years ago in a pioneering study A. T. Winfree [l| (see also 0) introduced a mean field model to describe 
T the limit-cycle behavior of large populations of biological oscillators. He discovered that systems of oscillators with 

''O ' randomly distributed frequencies remain incoherent when the variance of the frequencies is reduced, until a certain 
threshold is reached. Subsequently the oscillators begin to synchronize spontaneously and become locked. 
In its simplest form the model is defined by the set of equations (i = 1, ...iV, A^ ^ 1 ) 



k ^ 



^,(i)=c., + -^P(0,)i?(0,)- (1) 



3= 



m 

^^ Oiit) is the phase of the i—ih. oscillator; {tOi} describes a set of natural frequencies taken randomly from a distribution 

^p. \ g{uj). We shall assume below : g{u!) = I/27 for 76 [1 — 7, 1 + 7], gito) ~ otherwise; R{6i) is the sensitivity function 

• ' giving the response of the i— th oscillator; P{6j) is the influence function of the j—th oscillator. A common choice is 

O; R(e)^-sme, p(e)^i + cose . (2) 

00 , 

~. ' Despite its historical merits the Winfree's model has its own limitations. On one side it is complex enough not to 
^ . admit a full analytical treatment. On the other side it is not sufficiently sophisticated as to allow the treatment of 
realistic systems. Limitations of the former type were overcome by the work of Kuramoto [3| (for a review see[J|), who 
presented a model of oscillators related to the Winfree's model (it is the weak coupling limit of it) and analytically 
C^ ■ solvable in the mean field approximation. The Kuramoto's approach generated an intense theoretical work [5|, also 
motivated by the fact that the phenomenon of mutual synchronization of coupled nonlinear oscillators is ubiquitous in 
nature, with applications to neural networks, networks of cardiac pacemaker cells, populations of fireflies and crickets 
Q, [a| as well as arrays of Josephson junctions [3|- 

It must be also mentioned that, in spite of the complexity of the Winfree model, its phase diagram was object of 
investigation [s, ^ and in the (fc, 7) plane a rich structure was found. As to the versatility of the Winfree's model it 
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can be mentioned that it can describe different sets of pulse-coupled biological oscillators, see e.g. [lOi lllli H, 113 ■ I^ 
view of its relevance it might be interesting to look for extensions of the model that allow a different treatment of the 
couplings. Hopefully these extensions should enlarge the class of physical instances where the model can be usefully 
applied, for example experimental setups like arrays of Josephson's junctions [14 , Il5l | or the crowd synchronization 
phenomenon on the Millenium Bridge VE^. 

The purpose of this work is to generalize the Winfree model to the case of a global coupling depending on the fraction 
of synchronized oscillators. In a recent work [17| a modification of the original Kuramoto model was presented. The 
authors of 17] noted that the natural control parameter of the Kuramoto model is a coupling strength, analogous to 
the k parameter in ([T]), independent of the number of oscillators that are locked in frequencies. They suggested to 
generalize the theoretical model by allowing a dependence of the coupling on the number of locked oscillators. The 
technical way to achieve this is to introduce a functional dependence on the Kuramoto order parameter r, which is 
defined by the equation 



1 ^ 
^(i)e»'A(*)^ A^e'^'^W . (3) 

Clearly in a locked or partially locked phase r does not vanish and its variability produces the desired functional 
dependence. As the Kuramoto model is the averaged system of the Winfree model, it is desirable to study the effect 
in this latter case. Therefore in this paper we study the effect of an r-dependent coupling in the Winfree model, 
along the lines of [a, 0]. Our main result is that in some cases the modification produces an enlargement of the 
region, in the parameter space, where locking or partial locking of the oscillators is possible. It is worth noting that 
another interesting extension of the Kuramoto model, where additional powers of the order parameter are introduced 
to account for the dependence of the form of the coupling on its magnitude, has been recently studied [1^. The plan 
of the paper is as follows. In the next section we introduce the generalized Winfree's model, while in section 3 we 
describe how we obtain the transition lines in the phase diagrams. Some conclusions are drawn in section 4. 

II. THE GENERALIZED WINFREE'S MODEL 

We will consider the model defined by the set of equations 

^,(t)=c.. + ^^5]P(0,)i?(0O (4) 

i=i 

analogously to ([1]), with z a real parameter. It describes a set of coupled non linear oscillators, with coupling constant 
k r^~^. For z — I the model reduces to the Winfree model of H]). 

For N — > cx), the sum over all oscillators in ([T]) can be replaced by an integral, yielding the following equation for 
the velocity v — 6: 



where 



v{uj,0,t) = CO -a{t) sm0 (5) 



r27r /•1+7 

cr(i) = /cr^-M / {l + coae)p{e,t,Lj)g{uj)dwde . (6) 



Here p (w, 6, t) denotes the density of oscillators with phase 9 at time t. It satisfies the continuity equation 

dp d{pv) 

and the normalization condition 

I.2TT 



(7) 



dep{e,t,uj) = 1 . (8) 



for all Lu and any time. 

The phase diagram of the model with z — 1 was studied by Ariaratnam and Strogatz @]. They found a rich 
structure, comprising: i) locked phase (Lk), characterized by a common average frequency, ie a common value for the 



rotation number pi — \m\t^ooOi{t)/t] ii) partial locking (PL), characterized by macroscopic fractions of locked and 
unlocked oscillators; iii) incoherence (In), where no macroscopic fraction of oscillators is locked to a common frequency; 
iv) death (Dt) characterized by pi = for any i; v) partial death (PD), where only a fraction the pi vanishes. Moreover 
they found several hybrid states that can be seen as different realizations of the partial locking phase 1^ . 

We have numerically studied the model ([4]) with various values of the parameter z in the interval (0.5, 2.0) and 
various values of N , up to A^ = 1000, starting from a random initial configuration. In general stability in the results 
is found after 500 time steps. Our results are qualitatively similar to those of the original Winfree's model (z = 1), 
but the boundaries between the different phases depend on the actual value of z. Our numerical analysis is in general 
confirmed by the analytical study, see the next section. There is one case, the boundary locking / partial locking, 
when analytical results are not available, and the transition line must be evaluated only numerically. 

The main outcome of the study is that the value of the parameter z controls the size, in the phase diagram, of 
the incoherent phase: one may make the incoherent behavior less typical by choosing z < 1. The original Winfree's 
model, corresponding to z = 1, seems to be a special case. As we describe in the following section, we find that the 
partial locked region disappears for z > 1 and that, at fixed k and varying 7, the stability boundary of the locked 
phase corresponds to a continuous transition for z < 1 and first-order for z > 1. 

III. ANALYTICAL AND NUMERICAL RESULTS 

We shall define the transition lines between different regions in the phase diagram starting from areas where the 
solutions are stationary. Therefore we shall search for solutions characterized by a density po {^1^) f^d a velocity 
Wo {^tO) independent of time. Clearly also r and a are time independent in ([3]) and ([5]). 

The continuity equation ^ has stationary solutions; they satisfy po^o — C!{lo). From (O we see that \i lo < a one 
has the solution C{uj) = and therefore wq = 0. This implies that 

Pa = 5{e-e*), sinr = -, (9) 

a 

with the condition 

CT>l + 7. (10) 

The solution corresponds to the state of death (all the oscillators blocked at a fixed value of 0). We shall assume 
9* € (0, 7r/2), i.e. the result of the Winfree model [g| as we have numerically tested that this result holds also for 
generic z. We shall discuss this solution in subsection IIIIAI 
If ix) > cr, then C{uj) 7^ and we get 

Po(^,^) = ^^^. (11) 

w — cr sm 

From the normalization condition one has 



CH = ^^^ . (12) 

In the following we derive the stability boundaries between the phases of the model, generalizing the methods in @, 0] . 

A. Stability boundaries of the death phase 

From ([31) we get that, in general, 

rsinip = d9 dujgiuj) sin 6 po [uj ,9) 

Jq Jl-j 

/•2iT (•1+7 

rcosi/' = / d9 dujg{u!) cos 9 p{lu ,9) (13) 

Jo J1-7 



and in the death state 



rsinip — — (14) 



a 

f'l+7 



/•i+7 / ^2 

rcosip = / duj g{uj) \ 1 ^ (15) 



that can be employed to determine r and ij). 
Let us use ^ to get 



and adopt the definition 



kr 



— = 1 + J" cost/) , (16) 



G^(CT)=ari-^ . (17) 

An exphcit formula for the r.h.s of p^ is obtained by Eq. (|15p . One gets 



l+rcos'0 = H 

47 



i±^ jl - r^i^V - i^. 1 - ri^V + arcsin i±:i - arcsin i^ 



F,{a) . (18) 



The properties of F~f{a) were studied in [9]. For completeness we report here these results. It turns out that, as a 
function of a and for fixed 7 , F^[a) is a non-negative, increasing function, having the concavity down. From (|14p we 
also have r as a function of 7 and a: 



1 



^=\/^ + [^7(0-1] • (19) 



We distinguish the cases of large and small 7, the two ranges being separated by a limiting value 7^ £ (0, 1) that 
depends on z. For z = 1, 7^ = 0.2956 Q; let us generahze this result using a procedure similar to that of Q. At the 
same time we will characterize the two regions 7 < 7d and 7 > 7^. 

For a fixed value of 7 such that 7 < 7^, the two functions G^{a)/k and F^{a) can have one, two or no intersection, 
depending on the value of k. The value of a where the two curves are tangent, i.e. (Jd{l), satisfies (Jd{'j) > 1 + 7 and 
is the smallest value of a such that (|16p is satisfied. Therefore it characterizes the boundary. We can use P^ and 
the tangency condition 

^ = F», (20) 

to extract crd(7), getting rid of k. In this way one gets the boundary in the form 

The procedure can be repeated for various values of the parameter z characterizing the generalized Winfree's model. 
Increasing 7, Ud decreases and eventually it reaches the value 1 + 7. For any 7 such that 7 > 7d the values k are 
obtained by 



G^(l + 7) 



^7(1+7) 
The limiting value 7^ is the solution of the equation 

G-,(l+7) ^ G;(1 + 7) 
i^^(l + 7) F'AX^^) 



(7 > 7d) . (22) 



(23) 



The result, for various values of z in the interval (0, 2), is reported in FiglTJ 

The separation lines between the death region and the other phases (partial death, incoherence and locking- 
partial locking) are reported in Fig. [2] for four values ofz: 0.5, 1, 1.5 and 2. As stated above the analytical 
results are confirmed by the numerical analysis. It can be noted that the boundaries of the death phase are almost 
independent on z. This follows from the fact that the numerical values of r are quite close to unity, as can be seen 
expanding in the variable 7: 




0.25 



FIG. 1: The limiting value '^d as a function of z. The Winfree model is obtained for z = 1. 

B. Transition Incoherence/Partial Death 

Let us approach the boundary between these two phases from the the incoherence side. We use (|11|) in ^ . Since 
one must have cr < 1 — 7 the boundary is obtained putting a — 1 — 7. One has 



r cos "0 = 



1+7 



rsvaip = / du! g{u;)C{LL!) / d6 



1-7 



smf 



LU — (7 sm ( 



so that 



where 



r(7, a) = 



27 



f['-^)-f^'-' 



/(^) 



X \J x^ — 1 1 



2 2 

It follows that the boundary between the two regions is given by 

fc = (1 - 7) • r(7, 1-7) 



+ - In ix + \Jx^ - 1 ) 



(25) 



(26) 



(27) 



(28) 



with r given by ([26]) . Also these results are reported in Fig. [2l We note that increasing z the incoherent region 
increases while the partial death region decreases. 



C. Transition Incoherence/Partial Locking 

In order to determine the transition line we generalize the results of [9] to the case z 7^ 1. One adds to the static 
solution given in \\\\ : po {(^ ,0) = C{uj)/{lu — co sin 6), a small time-dependent perturbation: 



p{uj,d)^Po{oJ,9) + eTj{uj,9,t) . 



with e = 0^ and 






d0r]{LU,9,t) = 0. 



(29) 



(30) 
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FIG. 2: Phase diagram of the generahzed Winfree model. Results are for z — 0.5 (top-left), z = 1 (top-right), z = 1.5 (bottom- 
left) and z = 2 (bottom-right). Dt = death phase, PD = partial death, In = incoherence, Lk = locked phase, PL = partial 
locking. The phase diagram at z = 1 coincides with that depicted in [0|. In the case z = 1.5 the stability boundary between 
partial locking and incoherence has been drawn dashed-dotted because it is not observed in simulation (for z > 1 the partial 
locking region is absent, see the text). 



Similarly 



with 



a = (To + eai, 



(Tq = kr^ 



and 



0-1 



//'27r 
dujg{uj) / der]{u;,9,t), 



so that 

V = vq — eai {t) sin 9 , 
with Wo = a; — (To sm9. From the continuity equation one has, at first order in e: 



Searching solutions in the form 



dr] d{rivo) ujC{uj)cost 



77 = e^'M^,0), 



(31) 
(32) 

(33) 

(34) 

(35) 

(36) 



one finds 



d(hvo) u;C(uj)cosl 

Ah + — A t: 



86 



with 



A = kr' 



"^ / dhjg{uj) / 



deh{uj,0), 



wtiose solution is Q 



h(uj,9) ~ ^rla + bcos9 + csin^l, 



with 



a^ - A 



Therefore from (l38|) one gets 



with 



and, as in previous equation (j26p 



crowC'(w) 
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A2 + c^2 _ ^2 ' 
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= CTo/cro 
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= fcr^-1 







/ffo = / dujg{u;)ujX 



A2 + 6^2 _ ^2 
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f['-^)-f^'-' 



c^o 



CTO 



(37) 



(38) 



(39) 



(40) 



(41) 
(42) 



(43) 



(44) 



The transition line is obtained by taking the limit 3fi(A) — > 0^. Approximate results can be obtained performing a 
perturbative expansion for small 7, taking into account that in this limit also ctq ~^ 0- At the lowest order in 7 we 
can write A = i + Ai with Ai real. A better approximation can be obtained going up to the fifth order. The result is 

identical to the one found in [8|, [9| for z = 1, with the substitution k —^ ag. At the lowest order in 7 one gets r = — 
and therefore from (l42l) 



o-Q 



.1^- 



which shows that there are no solutions for z > 2. From (f43|) one finds 



87 ,, 

CTo = 1 1 

n 



I672 16{tt^ + 80)r 



and k is given by 



fc = 2 



m 



Oh') 



(45) 



(46) 



(47) 



In Fig. [2] the separation line between incoherence and partial locking is computed using the exact expression of 
Eqns. (|4T|) . (|42l) . and (|43|) . The approximate formula based on the expansion (|46|) is valid within 4 % for values of 
7 not larger than ~ 0.21 and z = 1; for z = 0.5 the validity is within 4 % for 7 < 0.19. It should be noted that for 
z > 1 moving from the right to the left at fixed k in the diagrams, one encounters the locking region before reaching 
the partial locking phase. Therefore, for z > 1 the partial locking region is basically absent, which means that, 
starting from a random initial configuration, the system never reaches a partial locking state. We remark that the 
occurrence of this phenomenon does not depends on the choice of ^(u;) uniform, indeed we verify that it holds also for 
a Lorentzian distribution. We refer the reader to |l3| for a discussion about the role of the shape of the distribution 
of frequencies in the Kuramoto model. 
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FIG. 3: The order parameter (time averaged) r versus 7, for low 7 and fixed k — 0.4, is reported for 2 — 0.5 (top-left, continuous 
transition between locked phase and partial locked phase) , z — 1 (top-right) z — 1.5 (bottom-left) and z — 2 (bottom-right). 
In the last three cases there is a first order transition between the locked phase and the incoherence phase. 

D. Transitions locking/partial Locking (z < 1) and locking/incoherence (z > 1) 

To derive the boundary from the locking to the partial locking phase, one should define the latter. This characteri- 
zation can be only heuristic, given the composite nature of the partial locked phase. We have obtained the transition 
line by the numerical solution of (|4]), using N — 800 oscillators and T = 1000 time units. To study the stability of 
the locked phase one should distinguish two cases, z < 1 and z > 1. In the case z < 1 we find a continuous transition 
from the locked phase towards the partial locking phase in which the time average of the order parameter decreases 
continuously between two limit values (see figure [3]). For z > 1 there is a first order transition from the locked phase 
to the incoherent phase as the time average of the order parameter r jumps from r ^ I (locked phase) to r < 0.1 
(incoherent phase), as shown in figure [31 In order to derive the curves limiting the locked phase, drawn in figUl we 
fix the coupling constant k; for each value we consider the plot of r versus 7 and (i) for 2 > 1 find the value of 7 at 
which the jump occurs (ii) for z < 1 find the value of 7 at which the curve has a flex point. The boundary line of 
the locked phase slightly depends on z (in the locked phase r « 1): as z decreases the locked phase region is slightly 
enlarged. 



IV. CONCLUSIONS 



We have presented a modification of the Winfree model to account for effective changes in the coupling constant 
among oscillators, as suggested by experiments on Josephson junctions and mechanical systems. The modification can 
be parametrized by a real number z. The case z — 1 corresponds to the Winfree model; z < 1 leads to a coupling which 
decreases as the order parameter r increases, thus enforcing the coupling at low r (low amount of synchronization); 



at z > 1 the coupling increases with r, i.e. the coupUng is weakened at low r. Using both an analytical approach and 
numerical simulations we have outlined the phase diagram of the model as z varies. 

As figure [2] clearly shows, the death phase region is almost independent on z, whilst the region of incoherence is 
strongly influenced by this parameter: for z < 1 it shrinks, as the effective coupling is strengthened at a low amount 
of synchronization, whereas it widens at z > 1 at the expenses of the partial death region. 

As far as the partial locking phase is concerned, we find that it disappears at z > 1 thus leading to the following 
phenomenon. At low k, as 7 is increased the system leaves the locking phase through a continuous transition (in the 
order parameter r) for z < 1, whilst for z > 1 the system undergoes a discontinuous transition while leaving the locked 
phase. This happens because for z < 1 the partial locking phase separates the locking and the incoherence phases, 
whereas for z > 1 the transition is directly onto the incoherence phase. The standard case z = 1, hence, appears to be 
rather special. It is worth noting that a similar change in the nature of the transition was noticed in the generalized 
Kuramoto model [17[ , and a discontinuous transition was experimentally seen in the synchronization of over-damped 
Josephson junctions [l5|, where physically the parameter z corresponds to the degree of feedback provided by a 
coupling resonator. Our results suggest strategies to control incoherent behavior in systems of interacting oscillators 
with coupling depending on the fraction of synchronized sub-units. 
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